rm(list=ls())
library(MASS)
library(stargazer)
library(knitr)
library(MASS)
library(foreign)
library(lme4)
library(car)

setwd("C:\\Users\\Daniel\\Dropbox\\Journal (II)\\Replication Files\\Karim\\")


fp <- read.csv("C:\\Users\\Daniel\\Dropbox\\Journal (II)\\Replication Files\\Karim\\PIWP_DataSet_2015 copy.csv")

##### NOTES on the Dataset
#The data was manually entered by enuemrators of the Center for Applied Research and Trainining in Liberia 
#The data were checked by the author periodically for errors.  If there were errors, the CART enumerators went back to the original surveys to correct the errors
#Some variables were duplicated and the 98 and 99 (I dont know and refuse to answer) were removed.  Thus, for some variables, there will be two columns one of the original responses and one with the I dont know and refuse to answer removed.  

#### Summary Statistics of peacekeeping contact 
summary(as.factor(fp$femaleunmilonly)) #100 contact only with female peacekeepers = 8% 
summary(as.factor(fp$maleunmilonly))  #144 contact with male peacekeeper = 11% 
summary(as.factor(fp$nounmilcontact))  #244 have had contact with peacekeepeer, 81% have not had contact with peacekeeper, 
#so 19% have had contact with peacekeeper 

#########################   RESULTS ########################

########Conact with Female Peacekeepers Table 1 Model 1 
########  Do you think that female peacekeepers are better than male peacekeepers?
#NOTE: morefemale.unmil.s is the variable without 98 or I dont know 

g<-fp$morefemale.unmil.s ~ fp$femaleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence    + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado  +fp$aflcontact
mod.logit1 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit1)

coeff.logit <- coefficients(mod.logit1)
provlnp233<-c(1,0,0,1,32.76328,0,0,1,0,0,0,0)
xb.provlnp233<-provlnp233%*%coeff.logit
predict.prob.win <- exp(xb.provlnp233)/(1+exp(xb.provlnp233))
predict.prob.win
# 0.1441589 (0.08402515-0.2381465)
provlnp233<-c(1,1,0,1,32.76328,0,0,1,0,0,0,0)
xb.provlnp233<-provlnp233%*%coeff.logit
predict.prob.win <- exp(xb.provlnp233)/(1+exp(xb.provlnp233))
predict.prob.win
# 0.2680937 (0.15364350-0.4305131)

mod.logit <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit)
mod.logit <- glm(formula=g,data=fp,family=binomial(link="logit"))
Z<-cbind(1,seq(0,1,0.01),0,1,32.76328,0,0,1,0,0,0,0)
Zb <- Z%*%coefficients(mod.logit) 
P <- exp(Zb)/(1+exp(Zb)) 
b.sim <- mvrnorm(1000,coefficients(mod.logit),vcov(mod.logit)) 
Zb.sim <- Z %*% t(b.sim) 
P.sim <- exp(Zb.sim)/(1+exp(Zb.sim))  
P.sim.graph <- t(apply(P.sim,1,sort)) 
P.sim.graph[,975]
P.sim.graph[,25]

########  Do you think that rape is a problem in Liberia?
# rapeproblem.2 is the variable without I dont know or 98 
g<-fp$rapeproblem.2 ~ fp$femaleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence    + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado +fp$aflcontact
mod.logit2 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit2)


######## Do you think that women should join the AFL?
######## Do you think that women should join the LNP?
#NOTE:Womenjoinafl.2 and womenjoin.lnp.2 is the variable without  I dont know or 98 

g<-fp$Womenjoinafl.2 ~ fp$femaleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence    + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado  +fp$aflcontact
mod.logit3 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit3)

g<-fp$womenjoin.lnp.2 ~ fp$femaleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence   + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado  +fp$aflcontact
mod.logit4 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit4)

######## Create Coefficient Table (Table 1 Model 1)
stargazer(mod.logit1,mod.logit2,mod.logit3,mod.logit4, type="html",  
          dep.var.labels=c("Fem. PK Better","Rape Problem","Women Join Mil.","Women Join Police"), 
          covariate.labels=c("Contact with Female Peacekeeper","No contact with peacekeeper ","Female","Age","Victim of Crime","Experienced Violence in War","Better off with UNMIL","Armed Group Member", "Contact with Female Police","Born in Monrovia","Contact with AFL") , 
          out="models.html", single.row = TRUE,digits = 2, title = "Table 1")


######## Subsets Sample Only: Table 2 and 3   
fem <- subset(fp, fp$Sex == "1" , select=c(hcode:whichdialect))
men <- subset(fp, fp$Sex == "0" , select=c(hcode:whichdialect))

with(fem, table(morefemale.unmil, femaleunmilonly))
with(fem, table(morefemale.unmil, maleunmilonly))
with(men, table(morefemale.unmil, femaleunmilonly))
with(men, table(morefemale.unmil, maleunmilonly))


#with(fem, table(morefemale.unmil.s, femaleunmilonly))
#with(fem, table(morefemale.unmil.s, maleunmilonly))
#with(fp, table(morefemale.unmil.s, femaleunmilonly))
#with(fp, table(morefemale.unmil.s, maleunmilonly))
#with(fp, table(Sex, femaleunmilonly))
#with(fp, table(morefemale.unmil, femaleunmilonly))
#with(fp, table(morefemale.unmil, maleunmilonly))


#######  Figure 1: Perceptions of Protection from Rape Provided by Different Actors

#Who do you think can best protect you from being raped? 
#rape2fp = prefer female peacekeeper to respond to rape 
#rape1mp = prefer male peacekeeper to respond to rape
#rape2femalepolice = prefer female officer to respond to rape 
#rape2malepolice = prefer male officer to respond to rape 
#rapet2femalecomm = prefer female community leader to respond to rape 
#rape2malecomleader = prefer male community leader to respond to rape 
summary(as.factor(fp$rape2fp))  #50/1280 = 0.04
summary(as.factor(fp$rape1mp))  #72/1280 = 0.06
summary(as.factor(fp$rape2femalepolice))  #505/1280= 0.39
summary(as.factor(fp$rape2malepolice))  #852/1280 = 0.67
summary(as.factor(fp$rapet2femalecomm)) #259/1280 = 0.20
summary(as.factor(fp$rape2malecomleader)) #534/1280 = 0.42


barplot(matrix(c(0.04,0.06,0.39,0.67,0.20,0.42),nr=2), ylab="Who do you think can best protect you from being raped?”  ", col=c("black","darkgray"),names.arg =c("Peacekeepers","Police","Community Leader"), ylim=c(0,0.8),beside=TRUE)
box()
legend("topleft", c("Female","Male"), pch=15,  col=c("black","darkgray"),  bty="n")


#########################   APPENDIX ########################


#Appendix Balance Table 

#Cognitive Ability: Balnaced 
summary(as.factor((fp$correctrace[fp$femaleunmilonly==1]))) #0.66
summary(as.factor(fp$correctrace[fp$femaleunmilonly==0]))  #0.4974446
summary(as.factor(fp$correctrace[fp$nounmilcontact==1]))  #0.4864078
summary(as.factor(fp$correctrace[fp$maleunmilonly==1])) #0.5763889
table <- matrix(c(66,83,34,61),ncol=2,byrow=TRUE)
colnames(table) <- c("Contact with Woman","Contact with Man")
rownames(table) <- c("Got it right","Didnt get it right")
table <- as.table(table)
table
prop.test(table, correct=FALSE)

#Age: Balnaced 
summary(fp$age[fp$femaleunmilonly==1]) # 35.34
summary(fp$age[fp$femaleunmilonly==0]) #32.51
summary(fp$age[fp$nounmilcontact==1]) #32.35
summary(fp$age[fp$maleunmilonly==1]) #33.67  
t.test(fp$age[fp$femaleunmilonly==1],fp$age[fp$maleunmilonly==1]) 

#Education: Balanced 
summary(fp$Education[fp$femaleunmilonly==1]) # 4.08
summary(fp$Education[fp$femaleunmilonly==0]) #3.5 
summary(fp$Education[fp$nounmilcontact==1]) #3.435 
summary(fp$Education[fp$maleunmilonly==1]) #3.972
t.test(fp$Education[fp$femaleunmilonly==1],fp$Education[fp$maleunmilonly==1]) 

#Income: Balanced 
summary(fp$Income[fp$femaleunmilonly==1]) # 4.83 
summary(fp$Income[fp$femaleunmilonly==0]) #5.059 
summary(fp$Income[fp$nounmilcontact==1]) #5.035
summary(fp$Income[fp$maleunmilonly==1]) #5.229 
t.test(fp$Income[fp$femaleunmilonly==1],fp$Income[fp$maleunmilonly==1]) 

#Household number: Unbalanced 
summary(fp$householdno[fp$femaleunmilonly==1]) #  4.44 
summary(fp$householdno[fp$femaleunmilonly==0]) #5.535 
summary(fp$householdno[fp$nounmilcontact==1]) #5.473
summary(fp$householdno[fp$maleunmilonly==1]) #5.979 
t.test(fp$householdno[fp$femaleunmilonly==1],fp$householdno[fp$maleunmilonly==1]) 

#Household number: Unbalanced 
summary(fp$childrennumber[fp$femaleunmilonly==1]) #   2.494
summary(fp$childrennumber[fp$femaleunmilonly==0]) #2.738 
summary(fp$childrennumber[fp$nounmilcontact==1]) #2.667
summary(fp$childrennumber[fp$maleunmilonly==1]) #3.248 
t.test(fp$childrennumber[fp$femaleunmilonly==1],fp$childrennumber[fp$maleunmilonly==1]) 

#Can Read: 
summary(fp$Read[fp$femaleunmilonly==1]) # 2.13 
summary(fp$Read[fp$femaleunmilonly==0]) #1.937 
summary(fp$Read[fp$nounmilcontact==1]) #1.914 
summary(fp$Read[fp$maleunmilonly==1]) #2.104 
t.test(fp$Read[fp$femaleunmilonly==1],fp$Read[fp$maleunmilonly==1]) 

#Sex: Not Balanced  
summary(as.factor((fp$sex[fp$femaleunmilonly==1]))) #0.15
summary(as.factor(fp$sex[fp$femaleunmilonly==0]))  #0.6062925
summary(as.factor(fp$sex[fp$nounmilcontact==1]))  #0.6302033
summary(as.factor(fp$sex[fp$maleunmilonly==1])) #0.4335664
table <- matrix(c(15,62,85,81),ncol=2,byrow=TRUE)
colnames(table) <- c("Contact with Woman","Contact with Man")
rownames(table) <- c("Female","Male")
table <- as.table(table)
table
prop.test(table, correct=FALSE)
# Men more likely to have contact with female peacekeeper 

#Victim Now: Unbalanced 
summary(as.factor((fp$nowvictim[fp$femaleunmilonly==1]))) #0.65
summary(as.factor(fp$nowvictim[fp$femaleunmilonly==0]))  #0.41
summary(as.factor(fp$nowvictim[fp$nounmilcontact==1]))  #0.40
summary(as.factor(fp$nowvictim[fp$maleunmilonly==1])) #0.46
table <- matrix(c(64,67,35,77),ncol=2,byrow=TRUE)
colnames(table) <- c("Contact with Woman","Contact with Man")
rownames(table) <- c("Victim Now","Not a Victim Now")
table <- as.table(table)
table
prop.test(table, correct=FALSE)

#Live in West Point: Balanced 
summary(as.factor((fp$community[fp$femaleunmilonly==1]))) #0.47
summary(as.factor(fp$community[fp$femaleunmilonly==0]))  #0.50
summary(as.factor(fp$community[fp$nounmilcontact==1]))  #0.51
summary(as.factor(fp$community[fp$maleunmilonly==1])) #0.42
table <- matrix(c(47,61,53,83),ncol=2,byrow=TRUE)
colnames(table) <- c("Contact with Woman","Contact with Man")
rownames(table) <- c("Live in WP","Live in PI")
table <- as.table(table)
table
prop.test(table, correct=FALSE)


#Experience War Violence: Unbalanced 
summary(as.factor((fp$warviolence [fp$femaleunmilonly==1]))) #0.38
summary(as.factor(fp$warviolence [fp$femaleunmilonly==0]))  #0.41
summary(as.factor(fp$warviolence [fp$nounmilcontact==1]))  #0.399
summary(as.factor(fp$warviolence [fp$maleunmilonly==1])) #0.51
table <- matrix(c(38,74,61,69),ncol=2,byrow=TRUE)
colnames(table) <- c("Contact with Woman","Contact with Man")
rownames(table) <- c("War Violence","No War Violence")
table <- as.table(table)
table
prop.test(table, correct=FALSE)
#Those with war violenc experience more likely to have contact with male 


#Contact with LNP: Balanced
summary(as.factor((fp$lnpcontact [fp$femaleunmilonly==1]))) #0.81
summary(as.factor(fp$lnpcontact [fp$femaleunmilonly==0]))  #0.6734521
summary(as.factor(fp$lnpcontact [fp$nounmilcontact==1]))  #0.6628019
summary(as.factor(fp$lnpcontact [fp$maleunmilonly==1])) #0.75
table <- matrix(c(81,108,19,36),ncol=2,byrow=TRUE)
colnames(table) <- c("Contact with Woman","Contact with Man")
rownames(table) <- c("Contact with LNP","No Contact with LNP")
table <- as.table(table)
table
prop.test(table, correct=FALSE)

#Contact with AFL: Unbalanced
summary(as.factor((fp$aflcontact [fp$femaleunmilonly==1]))) #0.84
summary(as.factor(fp$aflcontact [fp$femaleunmilonly==0]))  #0.3981324
summary(as.factor(fp$aflcontact [fp$nounmilcontact==1]))  #0.3800774
summary(as.factor(fp$aflcontact [fp$maleunmilonly==1])) #0.5277778
table <- matrix(c(84,76,15,68),ncol=2,byrow=TRUE)
colnames(table) <- c("Contact with Woman","Contact with Man")
rownames(table) <- c("Contact with AFL","No Contact with AFL")
table <- as.table(table)
table
prop.test(table, correct=FALSE)

#Contact times with UNMIL: Unbalanced  
summary(fp$contacttimesunmil[fp$femaleunmilonly==1]) #1.8
summary(fp$contacttimesunmil[fp$femaleunmilonly==0]) #0.167 
summary(fp$contacttimesunmil[fp$nounmilcontact==1]) #0.00388
summary(fp$contacttimesunmil[fp$maleunmilonly==1]) #1.333
t.test(fp$contacttimesunmil[fp$femaleunmilonly==1],fp$contacttimesunmil[fp$maleunmilonly==1]) 

#Contact times with LNP: Balanced   
summary(fp$contactimeslnp[fp$femaleunmilonly==1]) # 1.06 
summary(fp$contactimeslnp[fp$femaleunmilonly==0]) #1.014
summary(fp$contactimeslnp[fp$nounmilcontact==1]) #0.9923
summary(fp$contactimeslnp[fp$maleunmilonly==1]) # 1.174
t.test(fp$contactimeslnp[fp$femaleunmilonly==1],fp$contactimeslnp[fp$maleunmilonly==1]) 


#Armed Group: Unbalanced
summary(as.factor((fp$armedgroup [fp$femaleunmilonly==1]))) #0.22
summary(as.factor(fp$armedgroup [fp$femaleunmilonly==0]))  #0.03084833
summary(as.factor(fp$armedgroup [fp$nounmilcontact==1]))  #0.0311284
summary(as.factor(fp$armedgroup [fp$maleunmilonly==1])) #0.02962963
table <- matrix(c(22,4,78,135),ncol=2,byrow=TRUE)
colnames(table) <- c("Contact with Woman","Contact with Man")
rownames(table) <- c("Armed Group","Not in Armed Group")
table <- as.table(table)
table
prop.test(table, correct=FALSE)
#Female ex-combatants more likely to contact with female peaeckeeper 

#Contact with Female LNP:  Unbalanced 
summary(as.factor((fp$femalelnpcontact [fp$femaleunmilonly==1]))) #0.60
summary(as.factor(fp$femalelnpcontact [fp$femaleunmilonly==0]))  #0.3994911
summary(as.factor(fp$femalelnpcontact [fp$nounmilcontact==1]))  #0.3990338
summary(as.factor(fp$femalelnpcontact [fp$maleunmilonly==1])) #0.4027778
table <- matrix(c(60,58,40,86),ncol=2,byrow=TRUE)
colnames(table) <- c("Contact with Woman","Contact with Man")
rownames(table) <- c("Contact with Female LNP","No Contact with Female LNP")
table <- as.table(table)
table
prop.test(table, correct=FALSE)
#Those that have contact with female LNP are more likely to have contact with female peacekeeper


#Read:   
summary(as.factor((fp$read [fp$femaleunmilonly==1]))) #0.60
summary(as.factor(fp$femalelnpcontact [fp$femaleunmilonly==0]))  #0.3994911
summary(as.factor(fp$femalelnpcontact [fp$nounmilcontact==1]))  #0.3990338
summary(as.factor(fp$femalelnpcontact [fp$maleunmilonly==1])) #0.4027778
table <- matrix(c(60,58,40,86),ncol=2,byrow=TRUE)
colnames(table) <- c("Contact with Woman","Contact with Man")
rownames(table) <- c("Contact with Female LNP","No Contact with Female LNP")
table <- as.table(table)
table
prop.test(table, correct=FALSE)
#Those that have contact with female LNP are more likely to have contact with female peacekeeper


#Montserrado: 
summary(as.factor((fp$Montserrado[fp$femaleunmilonly==1]))) #0.43
summary(as.factor(fp$Montserrado [fp$femaleunmilonly==0]))  #0.3474576
summary(as.factor(fp$Montserrado [fp$nounmilcontact==1]))  #0.3503861
summary(as.factor(fp$Montserrado [fp$maleunmilonly==1])) #0.3263889
table <- matrix(c(43,47,57,97),ncol=2,byrow=TRUE)
colnames(table) <- c("Contact with Woman","Contact with Man")
rownames(table) <- c("Mont","Not from Mont")
table <- as.table(table)
table
prop.test(table, correct=FALSE)
#Female  more likely to be born in Monteserrado 

#Christian: 
summary(as.factor((fp$Christian[fp$femaleunmilonly==1]))) #0.88
summary(as.factor(fp$Christian [fp$femaleunmilonly==0]))  #0.8777589
summary(as.factor(fp$Christian [fp$nounmilcontact==1]))  #0.877176
summary(as.factor(fp$Christian [fp$maleunmilonly==1])) #0.8819444
table <- matrix(c(88,127,12,17),ncol=2,byrow=TRUE)
colnames(table) <- c("Contact with Woman","Contact with Man")
rownames(table) <- c("Christian","Not Christian")
table <- as.table(table)
table
prop.test(table, correct=FALSE)

######## Appendix Table 1: Contact with Male Peacekeepers 

g<-fp$morefemale.unmil.s ~ fp$maleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence    + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado  +fp$aflcontact
mod.logit1 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit1)

g<-fp$rapeproblem.2 ~ fp$maleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence    + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado +fp$aflcontact
mod.logit2 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit2)

g<-fp$Womenjoinafl.2 ~ fp$maleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence   + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado  +fp$aflcontact
mod.logit3 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit3)

g<-fp$womenjoin.lnp.2 ~ fp$maleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence  + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado  +fp$aflcontact
mod.logit4 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit4)

stargazer(mod.logit1,mod.logit2,mod.logit3,mod.logit4, type="html",  
          dep.var.labels=c("Fem. PK Better","Rape Problem","Women Join Mil.","Women Join Police"), 
          covariate.labels=c("Contact with Male Peacekeeper","No contact with peacekeeper ","Female","Age","Victim of Crime","Experienced Violence in War","Better off with UNMIL","Armed Group Member", "Contact with Female Police","Born in Monrovia","Contact with AFL") , 
          out="models.html", single.row = TRUE,digits = 2, title = "Appendix Table 1")

######## Appendix Table 2: Why Join the Police/Military  

summary(as.factor((fp$lnpjoin[fp$Sex==1]))) #0.5402985
summary(as.factor((fp$lnpjoin[fp$Sex==0]))) #0.3306452

summary(as.factor((fp$afljoin[fp$Sex==1]))) #0.3544118
summary(as.factor((fp$afljoin[fp$Sex==0]))) #0.4007491

summary(as.factor((fp$whyjoinafl[fp$Sex==1]))) 
summary(as.factor((fp$whyjoinafl[fp$Sex==0])))

summary(as.factor((fp$whyjoinlnp[fp$Sex==1]))) 
summary(as.factor((fp$whyjoinlnp[fp$Sex==0])))


########Appendix Table 3: "I don't knows?"  

g<-fp$morefemaleunmildk ~ fp$femaleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence   + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado  +fp$aflcontact + fp$Community
mod.logit1 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit1)

stargazer(mod.logit1, type="html",  
          dep.var.labels=c("Fem. PK Better","Wom. Join Mil."), 
          covariate.labels=c("Contact with Fem. Peacekeeper","No contact with peacekeeper ","Sex","Age","Victim of Crime","Experienced Violence in War","Better off with UNMIL","Armed Group Member", "Contact with Female Police","Born in Monrovia","Contact with AFL","West Point") , 
          out="models.html", single.row = TRUE,digits = 2, title = " Appendix Table 2: Responded I Don't Know")

######## Appendix Table 4: Women Only Subset Sample : Results 

fem <- subset(fp, fp$Sex == "1" , select=c(hcode:whichdialect))

g<-fem$morefemale.unmil.s ~ fem$femaleunmilonly + fem$nounmilcontact  +fem$Age + fem$nowvictim + fem$warviolence   + fem$better.2 +fem$armedgroup +fem$femalelnpcontact +fem$Montserrado +fem$aflcontact
mod.logit1 <- glm(formula=g,data=fem,family=binomial(link="logit"))
summary(mod.logit1)

g<-fem$rapeproblem.2 ~ fem$femaleunmilonly + fem$nounmilcontact  +fem$Age + fem$nowvictim + fem$warviolence    + fem$better.2 +fem$armedgroup +fem$femalelnpcontact +fem$Montserrado +fem$aflcontact
mod.logit2 <- glm(formula=g,data=fem,family=binomial(link="logit"))
summary(mod.logit2)

g<-fem$Womenjoinafl.2 ~ fem$femaleunmilonly + fem$nounmilcontact  +fem$Age + fem$nowvictim + fem$warviolence   + fem$better.2 +fem$armedgroup +fem$femalelnpcontact +fem$Montserrado  +fem$aflcontact
mod.logit3 <- glm(formula=g,data=fem,family=binomial(link="logit"))
summary(mod.logit3)

g<-fem$womenjoin.lnp.2 ~ fem$femaleunmilonly + fem$nounmilcontact  +fem$Age + fem$nowvictim + fem$warviolence    + fem$better.2 +fem$armedgroup +fem$femalelnpcontact +fem$Montserrado  +fem$aflcontact
mod.logit4 <- glm(formula=g,data=fem,family=binomial(link="logit"))
summary(mod.logit4)

stargazer(mod.logit3,mod.logit4, type="html",  
          dep.var.labels=c("Women Join Mil.","Women Join Police"), 
          covariate.labels=c("Contact with Fem. Peacekeeper","No contact with peacekeeper ","Age","Victim of Crime","Experienced Violence in War","Better off with UNMIL","Armed Group Member", "Contact with Female Police","Born in Monrovia","Contact with AFL") , 
          out="models.html", single.row = TRUE,digits = 2, title = "Table 2: Women Only Sample")


####### Appendix 5: Community Level Fixed Effects 

g<-fp$morefemale.unmil.s ~ fp$femaleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence    + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado  +fp$aflcontact + fp$Community
mod.logit1 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit1)

g<-fp$rapeproblem.2 ~ fp$femaleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence    + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado  +fp$aflcontact + fp$Community
mod.logit2 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit2)

g<-fp$Womenjoinafl.2 ~ fp$femaleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence    + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado  +fp$aflcontact+ fp$Community
mod.logit3 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit3)

g<-fp$womenjoin.lnp.2 ~ fp$femaleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence   + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado  +fp$aflcontact + fp$Community
mod.logit4 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit4)


stargazer(mod.logit1,mod.logit2,mod.logit3,mod.logit4, type="html",  
          dep.var.labels=c("Fem. PK Better","Fem Respond Rape","Women Join Mil.","Women Join Police"), 
          covariate.labels=c("Contact with Female Peacekeeper","No contact with peacekeeper ","Female","Age","Victim of Crime","Experienced Violence in War","Better off with UNMIL","Armed Group Member", "Contact with Female Police","Born in Monrovia","Contact with AFL","West Point") , 
          out="models.html", single.row = TRUE,digits = 2, title = "Table 1")

####### Appendix 6: West Point and Peace Island Sample Only 

split <- subset(fp, fp$Community == "0" , select=c(hcode:whichdialect))

g<-split$morefemale.unmil.s ~ split$femaleunmilonly + split$nounmilcontact + split$Sex +split$Age + split$nowvictim + split$warviolence    + split$better.2 +split$armedgroup +split$femalelnpcontact +split$Montserrado  +split$aflcontact 
mod.logit1 <- glm(formula=g,data=split,family=binomial(link="logit"))
summary(mod.logit1)

with(split, table(morefemale.unmil.s, femaleunmilonly))
with(split, table(morefemale.unmil.s, maleunmilonly))

g<-split$rapeproblem.2 ~ split$femaleunmilonly + split$nounmilcontact + split$Sex +split$Age + split$nowvictim + split$warviolence    + split$better.2 +split$armedgroup +split$femalelnpcontact +split$Montserrado  +split$aflcontact 
mod.logit2 <- glm(formula=g,data=split,family=binomial(link="logit"))
summary(mod.logit2)

with(split, table(rapeproblem.2, femaleunmilonly))
with(split, table(rapeproblem.2, maleunmilonly))

g<-split$Womenjoinafl.2 ~ split$femaleunmilonly + split$nounmilcontact + split$Sex +split$Age + split$nowvictim + split$warviolence    + split$better.2 +split$armedgroup +split$femalelnpcontact +split$Montserrado  +split$aflcontact
mod.logit3 <- glm(formula=g,data=split,family=binomial(link="logit"))
summary(mod.logit3)

g<-split$womenjoin.lnp.2 ~ split$femaleunmilonly + split$nounmilcontact + split$Sex +split$Age + split$nowvictim + split$warviolence   + split$better.2 +split$armedgroup +split$femalelnpcontact +split$Montserrado  +split$aflcontact 
mod.logit4 <- glm(formula=g,data=split,family=binomial(link="logit"))
summary(mod.logit4)

stargazer(mod.logit1,mod.logit2,mod.logit3,mod.logit4, type="html",  
          dep.var.labels=c("Fem. PK Better","Fem Respond Rape","Women Join Mil.","Women Join Police"), 
          covariate.labels=c("Contact with Female Peacekeeper","No contact with peacekeeper ","Female","Age","Victim of Crime","Experienced Violence in War","Better off with UNMIL","Armed Group Member", "Contact with Female Police","Born in Monrovia","Contact with AFL") , 
          out="models.html", single.row = TRUE,digits = 2, title = "Table 1")


######### Appendix 8:  Male Only Sample 

fem <- subset(fp, fp$Sex == "0" , select=c(hcode:whichdialect))

g<-fem$morefemale.unmil.s ~ fem$femaleunmilonly + fem$nounmilcontact  +fem$Age + fem$nowvictim + fem$warviolence   + fem$better.2 +fem$armedgroup +fem$femalelnpcontact  +fem$Montserrado  +fem$aflcontact
mod.logit1 <- glm(formula=g,data=fem,family=binomial(link="logit"))
summary(mod.logit1)

mod.logit2 <- glm(fem$womenrape ~ fem$femaleunmilonly + fem$nounmilcontact  +fem$Age + fem$nowvictim + fem$warviolence   + fem$better.2 +fem$armedgroup +fem$femalelnpcontact +fem$Montserrado  +fem$aflcontact , family = "poisson", data = fem)
summary(mod.logit2)

g<-fem$Womenjoinafl.2 ~ fem$femaleunmilonly + fem$nounmilcontact  +fem$Age + fem$nowvictim + fem$warviolence   + fem$better.2 +fem$armedgroup +fem$femalelnpcontact +fem$Montserrado  +fem$aflcontact
mod.logit3 <- glm(formula=g,data=fem,family=binomial(link="logit"))
summary(mod.logit3)

g<-fem$womenjoin.lnp.2 ~ fem$femaleunmilonly + fem$nounmilcontact  +fem$Age + fem$nowvictim + fem$warviolence    + fem$better.2 +fem$armedgroup +fem$femalelnpcontact +fem$Montserrado  +fem$aflcontact
mod.logit4 <- glm(formula=g,data=fem,family=binomial(link="logit"))
summary(mod.logit4)


stargazer(mod.logit1,mod.logit2,mod.logit3,mod.logit4, type="html",  
          dep.var.labels=c("Fem. PK Better","Fem.PK Rape","Women Join Mil.","Women Join Police"), 
          covariate.labels=c("Contact with Fem. Peacekeeper","No contact with peacekeeper ","Age","Victim of Crime","Experienced Violence in War","Better off with UNMIL","Armed Group Member", "Contact with Female Police","Born in Monrovia","Contact with AFL") , 
          out="models.html", single.row = TRUE,digits = 2, title = "Table 3: Men Only Sample")



#############  BEBER ET AL SURVEY RESULTS 
rm(list=ls())
library(MASS)
library(stargazer)
library(knitr)
library(MASS)
library(foreign)
library(lme4)
library(car)

setwd("~/Desktop")

mon <- read.csv("~/Dropbox/II Solo/Final/Replication/Monrovia2012.csv")

#setwd("~/Dropbox/1325 project/Data")

#mon<-as.data.frame(read.dta("Data all.dta"),use.value.labels = TRUE)
#setwd("~/Desktop")
#write.csv(mon, file = "Monrovia2012.csv")

### CODEBOOK: 
#q19: WHEN was the last time you socialized with UNMIL PKO personnel?
#q21: Was this a male or a female peacekeeper?
#q25: Presence of UNMIL personnel has increased OWN personal security?
#savings = Do you have any savings?
#cognitive (q10) = You overtake the person in 2nd place in a race. Your position? Correct Answer 
#muslim = using the variable "religion (q5)" a dummy bariable was created based on whether the respondent said they were Muslim 
#q1age = age 
#sex = q2 or sex (Female = 1)
#wartrauma (q148): Were any years since start of war really hard times for you?
#psu: community fixed effects 

summary(as.factor(mon$psu))

summary(as.factor(mon$q19))  #21 per cent had contact in the last year 
summary(as.factor(mon$q21)) #  0.848 with a male and  0.152 with female 
summary(as.factor(mon$q25))

#Three new columns created using q21 of contact with female peacekeepers (q21_women), contact with male peacekeepers (mon$q21_men), and no contact with peacekeepers (mon$q21_none).
summary(as.factor(mon$q21_women))
summary(as.factor(mon$q21_men))
summary(as.factor(mon$q21_none))

g<-mon$q25 ~ mon$q21_women + mon$q21_none + mon$q1_age +  + mon$sex + mon$savings  + mon$muslim + mon$cognitive + mon$wartrauma + as.factor(mon$psu)
mod.logit1 <- glm(formula=g,data=mon,family=binomial(link="logit"))
summary(mod.logit1)

g<-mon$q25 ~ mon$q21_men + mon$q21_none + mon$q1_age +   mon$sex + mon$savings + mon$muslim + mon$cognitive + mon$wartrauma + as.factor(mon$psu)
mod.logit2 <- glm(formula=g,data=mon,family=binomial(link="logit"))
summary(mod.logit2)

g<-mon$q25 ~  mon$q21_women + mon$q21_men + mon$q1_age  + mon$sex + mon$savings +  mon$muslim + mon$cognitive + mon$wartrauma + as.factor(mon$psu)
mod.logit3 <- glm(formula=g,data=mon,family=binomial(link="logit"))
summary(mod.logit3)

stargazer(mod.logit1, mod.logit2, mod.logit3, type="html",  
          dep.var.labels=c("Has the presence of UNMIL increased OWN personal security?"), 
          covariate.labels=c("Contact with Female Peacekeeper","Contact with Male Peacekeeper","No contact with peacekeeper ","Age","Female","Savings","Muslim","High Cognitive Ability","War Trauma") , 
          out="models.html", single.row = TRUE,digits = 2, title = "Table 4: Representative Sample of Monrovia")

########################################################################################################################
########################################################################################################################
# Below Analysis is not included in the paper
##########Subset Sample by Sex 

fem <- subset(mon, mon$r_sex1 == "1" , select=c(id_hh:TIMESTAMPS))

g<-fem$q25 ~ fem$q21_women + fem$q21_none + fem$r_age + fem$r_working  + fem$no_contact_with_foreigner+ mon$muslim + fem$cognitive
mod.logit1 <- glm(formula=g,data=fem,family=binomial(link="logit"))
summary(mod.logit1)

g<-fem$q25 ~ fem$q21_men + fem$q21_none + fem$r_age + fem$r_working + fem$no_contact_with_foreigner+ mon$muslim + fem$cognitive
mod.logit2 <- glm(formula=g,data=fem,family=binomial(link="logit"))
summary(mod.logit2)

g<-fem$q25 ~  fem$q21_women +fem$q21_men + fem$r_age + fem$r_working  + fem$no_contact_with_foreigner+ mon$muslim + fem$cognitive
mod.logit3 <- glm(formula=g,data=fem,family=binomial(link="logit"))
summary(mod.logit3)

stargazer(mod.logit1,mod.logit2,mod.logit3, type="html",  
          dep.var.labels=c("Has the presence of UNMIL increased OWN personal security?"), 
          covariate.labels=c("Contact with Female Peacekeeper","Contact with Male Peacekeeper","No contact with peacekeeper ","Age","Have a Job","No Contact with Foriegners","Christian","Muslim","High Cognitive Ability") , 
          out="models.html", single.row = TRUE,digits = 2, title = "Table 4: Representative Sample of Monrovia (Women Only)")


fem <- subset(mon, mon$r_sex1 == "1" , select=c(id_hh:TIMESTAMPS))

g<-fem$q25 ~ fem$q21_women + fem$q21_none + fem$r_age + fem$r_working  + fem$no_contact_with_foreigner+ mon$muslim + fem$cognitive
mod.logit1 <- glm(formula=g,data=fem,family=binomial(link="logit"))
summary(mod.logit1)

coeff.logit <- coefficients(mod.logit1)
provlnp233<-c(1,0,0,35.33202,0,1,1,0)
xb.provlnp233<-provlnp233%*%coeff.logit
predict.prob.win <- exp(xb.provlnp233)/(1+exp(xb.provlnp233))
predict.prob.win
# 0.9592222
provlnp233<-c(1,1,0,35.33202,0,1,1,0)
xb.provlnp233<-provlnp233%*%coeff.logit
predict.prob.win <- exp(xb.provlnp233)/(1+exp(xb.provlnp233))
predict.prob.win
# 0.990872


g<-fem$q25 ~ fem$q21_men + fem$q21_none + fem$r_age + fem$r_working + fem$no_contact_with_foreigner+ mon$muslim + fem$cognitive
mod.logit2 <- glm(formula=g,data=fem,family=binomial(link="logit"))
summary(mod.logit2)

g<-fem$q25 ~  fem$q21_women +fem$q21_men + fem$r_age + fem$r_working  + fem$no_contact_with_foreigner+ mon$muslim + fem$cognitive
mod.logit3 <- glm(formula=g,data=fem,family=binomial(link="logit"))
summary(mod.logit3)

stargazer(mod.logit1,mod.logit2,mod.logit3, type="html",  
          dep.var.labels=c("Has the presence of UNMIL increased OWN personal security?"), 
          covariate.labels=c("Contact with Female Peacekeeper","Contact with Male Peacekeeper","No contact with peacekeeper ","Age","Have a Job","No Contact with Foriegners","Christian","Muslim","High Cognitive Ability") , 
          out="models.html", single.row = TRUE,digits = 2, title = "Table 5: Representative Sample of Monrovia (Men Only)")


#summary(as.factor((mon$q25[mon$q21_women==1]))) #94% think increased own security 
#summary(as.factor((mon$q25[mon$q21_women==0]))) #73% think increased own security 

#summary(as.factor((mon$q25[mon$q21_men==1]))) #81% think increased own security 
#summary(as.factor((mon$q25[mon$q21_men==0])))  #71% think increasen own security  

#test<-cbind(c(94),c(81))
#whenborn<-as.matrix(test)
#barplot(whenborn, beside=TRUE, main="Percentage of respondents that think UNMIL has increased personal security
        #(Monrovia)",names.arg =c("Contact with Female PKO","Contact with Male PKO"),ylim=c(0,100), col="dark green")



